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INTRODUCTION 


Predictions  of  gas  turbine  combustor  performance  and  pollutant  emission 
characteristics  require  modeling  procedxires  possessing  a high  degree  of  sophis- 
tication. ftLst  attempts  at  modeling  combustion  systems  have  been  largely 
frustrated  by  the  complexity  of  the  coupled  hydrodynamic  and  chemical  pro- 
cesses. The  difficulty  can  be  largely  attributed  to  the  lack  of  understanding 
of  the  flow  processes  which,  through  the  exchange  of  heat,  mass  and  momentvun, 
can  directly  relate  to  pollutant  formation  and  combxistion  efficiency.  In  ad- 
dition, a reliable  calculation  procedure  must  be  capable  of  properly  treating 
the  reverse  flow  recirculation  regions  occurring  in  combustion  devices.  The 
objective  of  the  present  investigation  was  to  extend  an  existing  and  relative- 
ly efficient  Uhited  Technologies  Research  Center  three-dimensional  Navier- 
Stokes  calculation  procedure  (Refs.  1 and  2)  so  that  it  would  be  able  to  com- 
pute combusting  flows.  In  particular,  the  procedure  developed  solves  the  time- 
averaged  three  dimensional  canpressible  Navier-Stokes  equations  including  a two- 
equation  transport  turbvilence  model,  a pseudo-kinetic  hydrocarbon  chemistry 
model,  a liquid  droplet  vaporization  model,  a single  frequency  radiation  model 
and  a finite  rate  nitrogen  oxide  chemistry  model.  The  intent  in  the  develop- 
ment of  the  flow  models  described  above  was  to  eliminate  luidue  complexity  and 
sophistication,  simultaneously  providing  a reasonably  good  framework  within 
which  refinements  could  be  easily  implemented  at  a futiore  time,  if  warranted 
by  comparisons  with  experimental  data. 

The  theoretical  analysis  pertaining  to  the  present  study  is  described  in 
detail  in  Volumes  I and  II  of  the  present  report.  In  view  of  the  preliminary 
nature  of  the  present  work,  this  user's  manual  is  not  intended  to  be  completely 
comprehensive.  A general  description  of  the  main  feat\ares  of  the  ccmputer 
program  is  provided  axong  with  a description  of  the  input  variables  and  a 
summary  of  the  printed  output. 


1.  Briley,  W.  R.  and  H.  McDonald:  An  Implicit  Numerical  Method  for  the 

Multidimensional  Compressible  Navier-Stokes  Equations,  United  Aircreift 
Research  Laboratories,  Report  M9H363-6,  November  1973. 

2,  Briley,  W.  R. , H.  McDonald,  and  H.  J.  Gibeling:  Solution  of  the  Multi- 

dimensional Compressible  Navier-Stokes  Equations  by  a Generalized  Implicit 
Method.  United  Technologies  Research  Center,  Report  R75-9II363-I5, 

January  1976. 
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COMPUTER  PROGRAM 


^tode  of  Operation 


The  computer  program  is  written  to  solve  the  three-dimensional  time- 
averaged  compressible  Navier-Stokes  equations  including' a two-equation  trans- 
port turbulence  model,  a psexido-kinetic  hydrocarbon  chemistry  model,  a liquid 
droplet  vaporization  model,  and  a single  frequency  radiation  transport  model. 
Also,  a finite  rate  nitrogen  oxide  chemistry  model  is  available  using  a restart 
procedure  to  be  described.  The  numerical  procedure  solves  the  time- dependent 
governing  equations  beginning  with  a specified  starting  solution  and  appropri- 
ate bovindary  conditions  until  a steady  state  solution  is  achieved. 

Several  coordinate  system  options  are  available  in  the  present  computer 
program.  Cartesian  coordinates  are  available  for  flow  problems  in  rectangular 
ducts  with  injection  ports  on  any  wall.  A generalized  axisymmetric  orthogonal 
curvilineeir  coordinate  system  is  available  for  flow  in  axisymmetric  combustors 
with  a discrete  circumferential  distribution  of  injectidn  ports.  However,  the 
computer  program  has  been  tested  only  for  cylindrical-polar  and  cartesian  co- 
ordinates, and  an  earlier  version  of  the  deck  (Refs.  1 and  2)  was  thoroughly 
evaluated  with  cartesian  coordinates.  The  necessary  link  to  generate  the  metric 
scale  factors  for  the  generalized  coordinate  system  has  not  been  established 
at  this  time . 

The  input  data  for  this  program  is  all  in  NAMELIST  format.  Reference 
values  are  required  for  velocity,  length,  density,  temperature,  specific  heat, 
dynamic  viscosity  and  molecular  weight.  Default  values  for  all  of  these 
quantities  except  velocity  and  length  are  assumed  in  the  program,  and  corre- 
spond to  air  at  standard  conditions.  In  addition,  a nomln- i nondimensional 
pressure  is  required  for  computation  of  a pressure  difference,  which  is  used 
in  the  momentum  conservation  equations  to  eliminate  truncation  errors  in  pres- 
sure gradient  terms.  The  specific  heat  ratio  (y)  is  needed  to  con^jute  a 
reference  fl*ozen  sound  speed  used  in  calculating  e frozen  flow  Mach  nvunber  for 
output  only.  The  units,  defaiilt  values,  and  Fortran  symbols  for  the  reference 
quantities  are  given  in  Table  I.  Default  values  for  other  input  variables 
are  given  in  the  input  description  subsection. 

The  actual  equations  to  be  solved  and  the  appropriate  boundary  conditions 
must  be  specified  via  input  cards  to  override  default  values.  When  a function 
boundary  condition  is  employed  at  a wall  or  port,  the  desired  boundary  value 
of  the  dependent  variable  must  be  specified  via  the  $BCDATA  namelist  input 
unless  it  is  identically  zero. 


The  initial  flow  field  estimate  may  be  obtained  in  several  ways.  First, 
the  initial  flow  field  can  be  generated  internally  using  only  inlet  profiles 
and  a specified  pressure  drop  from  input.  Alternatively,  the  entire  starting 
flow  field  can  be  specified  via  input,  if  sufficiently  smooth  consistent  flow 
field  data  is  available.  Finally,  the  initial  flow  field  may  be  read  from  a 
magnetic  tape  or  a disk  file  generated  by  a previous  run.  The  latter  mode  is 
used  for  restarting  an  ongoing  solution,  and  its  iise  is  highly  recommended  to 
allow  checking  the  solution  for  errors.  At  present  the  restart  procedure  is 
required  when  initiating  the  kinetic  nitrogen  oxide  chemistry  calculation  from 
a converged  reacting  flow  field  solution. 

Once  the  initial  and  botindary  data  have  been  specified,  the  program  begins 
computation.  The  calculation  will  proceed  for  either  a specified  number  of 
time  steps  or  until  the  solution  changes  by  less  than  a liser  specified  toler- 
ance between  successive  time  steps,  which  ever  occurs  first. 


The  Calculation  Procedure 

The  MINT  (Multidimensional  Implicit  Nonlinear  Time-dependent)  Combustor 
Code  ccmbines  a main  program  and  a series  of  subroutines  to  perfomn  the 
required  calculations.  A block  diagram  of  the  program  is  shown  in  Fig.  1. 

The  main  program  DAL  contains  comment  and  data  statements  pertaining  to 
the  use  of  the  program.  Also,  DAL  calls  the  main  sequence  of  overlayed 
subroutines.  The  input  is  read,  initial  parameters  are  set,  and  an 
initial  flow  field  is  calculated  in  subroutine  READA  and  the  associated 
subroutines  it  calls.  If  a nitrogen  oxide  chemistry  calculation  is  being 
started,  subroutine  QNOXI  is  called  to  initialize  all  variables  necessary  for 
that  calculation.  Then  the  execution  control  stibroutine  EXEC  is  called  from 
DAL.  Subroutine  EXTBV  is  called  from  EXEC  to  initialize  the  remaining 
quantities  needed  to  start  the  calculation,  and  the  initial  flow  field  is 
printed  out  by  PRNTA.  The  time  step  loop  is  then  entered  in  EXEC  and  sub- 
routine TSTEP  advances  the  solution  through  one  time  step.  Various  output 
functions  are  then  performed  in  EXEC  and  the  next  time  step  is  computed.  The 
solution  veuriables  are  printed  out  at  a specified  time  step  frequency  by 
subroutine  PRNTB.  The  procedure  then  has  completed  one  time-step  calculation. 
After  completion  of  the  specified  number  of  time  steps  or  convergence  to  a 
steady  solution,  the  control  is  ret\irned  to  the  main  program  DAL  where  sub- 
routine RESTRT  is  called  to  produce  a restart  dump  if  desired,  and  finally 
subroutine  PRNTF  is  called  to  obtain  the  final  variables  and  a limited  number 
of  printer  plots. 


General  Featvires  of  the  Program 


L 


The  program  is  written  in  FORTRAN  IV  and  may  require  one  or  more  disk 
files  and/or  tapes  for  temporary'  snd  permanent  data  storage. 


The  fluid  is  a compressible  gas  wi-th  variable 
thermodynamic  properties. 

Three-dimensional  rectangular  ducts  or  axis:,’mmetric 
geometries  with  discrete  injection  ports  ma.y  be 
considered  (Fig.  2). 

The  flow  may  be  turbulent  or  laminar  with  heat  and 
mass  transfer.  I^ydrocarbon  chemistry,  t'wo-phase 
flows,  and  flows  with  recirculation  may  be  treated. 
Effective  exchange  coefficients  and  a two-equation 
model  are  employed  to  model  the  t'urb’ulence.  A 
pseudo-kinetic  chemistry  and  a droplet  model  are 
used  for  the  chemical  and  two-phase  flow  modeling. 

Types  of  Initial  Flow  The  initial  flow  field  may  be  read  entirely  from 

Field  Estimate  cards  or  generated  internally  by  the  program  (cold 

start),  or  be  obtained  from  a magnetic  tape 
or  disk  file. 

Types  of  Boundary  Function  or  derivative  boundary  conditions  may  be 

Conditions  specified  for  any  dependent  variable.  The  time- 

dependent  continuity  equation  is  provided  as  the 
boundary  condition  for  the  continuity  equation  at 
a solid  wall.  At  an  inlet  port,  it  is  recommended 
that  the  density  be  specified  as  the  contin\iity 
equation  boundary  condition  and  a second- derivative 
boundary  condition  on  the  pressure  should  then  be 
employed  as  the  energy  equation  boundarj'  condition. 
In  any  case,  the  exit  port  boundary  condition  for 
the  continuity  equation  should  be  a specified  static 
pressure. 


Type  of  Fluid 


Type  of  Geometry 


Types  of  Flow 


Types  of  Printed  Output  An  abbreviated  summary  of  the  solution  is  printed 

out  after  each  time  step.  In  addition,  the  user 
may  select  a time  step  interval  at  which  a full 
print  of  the  flow  field  variables  is  printed.  A 
full  print  is  also  executed  after  the  final  time 


U 
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step  of  a rm  along  with  printer  generated  time 
history  and  spatial  plots  of  selected  flow  variables. 

The  IRERJN  Option  Using  this  option  the  user  may  have  a restart 

written  on  a magnetic  tape  or  a permanent  disk 
file.  This  file  contains  sufficient  informa- 
tion to  start  a new  run  commencing  where  the 
previous  one  ended.  The  restart  file  is  written 
after  the  calculation  has  proceeded  for  the  user 
specified  a number  of  time  steps,  or  at  a user 
specified  time  step  interval,  or  if  the  program 
is  in  danger  of  exceeding  the  computer  run  time 
limit. 


Namelist  Input  Description 


Namelist  or 
variable  name 


$READ1 


IREST 


I 


I 


IREPUN 


I0TAPE 

INTAPE 

I0TAP1 

INTAPl 


Description 

Restart  option  selection 


IREST  =0.  A new  calculation  is  being  started. 
Read  input  from  ciirds.  Default. 

IREST  = 1.  Restart  option.  Read  namelist  data 
from  l/O  \init  number  INTAPl  and  binary  (unfor- 
matted) array  data  from  unit  number  INTAPE. 

IREST  = U.  Restart  option.  Read  namelist  data 
from  unit  number  INTAPl  and  formatted  array 
data  from  unit  number  INTAPE.  In  this  case, 
INTAPl  may  equal  INTAPE. 

IREPUN  = 0.  No  restart  dtunp  is  written  during 
current  r\in. 

IREPUN  = 1.  Write  restart  dump  at  time  step 
interval  NTREST  (see  $READU  namelist).  Name- 
list  data  is  written  onto  unit  number  I0TAP1. 
Binary  (unformatted)  array  data  is  written  onto 
unit  n’umber  I0TAPE. 

IREHJN  = 4.  Same  as  IREPUN  = 1 option  except 
formatted  array  data  is  written  onto  unit  number 
I0TAPE,  In  this  case  I0TAP1  may  equal  I0TAPE. 

Output  unit  number  for  dependent  variable  array 
restart  data.  Default  is  I0TAPE  = 10. 

Input  unit  number  for  dependent  variable  array 
restart  data.  Default  is  INTAPE  = 9. 

Output  \uiit  number  for  namelist  restart  data. 
Default  is  I0TAP1  = 20. 

Input  unit  n\imber  for  namelist  restart  data. 
Default  is  INTAPl  = 19 . 
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Namelist  or 
variable  name 

NFINP 

NF0UT 

$READ2 

NUMDX 

NUMDY 

NUMDZ 

XMIN  or  XGMIN(l) 
XMAX  or  XGMAX(l) 
YMIN  or  XGMIN(2) 
YMAX  or  XGMAX(2) 


Description 
Not  used. 

Not  used. 

Computational  grid  selection 

Number  of  interior  grid  points  in  the  X - dir- 
ection. Two  additional  boundary  points^are  also 
required . 


Number  of  interior  grid  points  in  the  X - dir- 
ection. Two  additional  boundary  points^are  also 
required. 


Number  of  interior  grid  points  in  the  X - dir- 
ection. Two  additional  boundary  points'^are  also 
required . 

Nondimens ional  value  of  X at  boundary  number  1 
(Fig.  3);  identified  by  iAdex  NXl  ■ 1. 

Nondimens ional  value  of  X^  at  boundary  number  2 
(Fig.  3);  identified  by  index  NX2.  XMAX  > XMIN. 

Nondimens ional  value  of  X„  at  boundary  number  3 
(Fig.  3);  identified  by  index  NYl  ■ 1. 

Nondimens ional  value  of  X^  at  boundary  number  4 
(Pig.  3);  identified  by  index  NY2.  YMAX  > YMIN. 


ZMIN  or  XC31IN(3) 


ZMAX  or  XGMAX(3) 


Nondimensional  vsilue  of  X at  boundary  number  ^ 
(Fig.  4);  identified  by  index  NZl  ■ 1. 

Nondimensional  value  of  X^  at  boundary  number  6 
(Fig.  4);  identified  by  index  NZ2.  ZMAX  > ZMIN. 
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Namelist  or 
variable  name 


IGRID(I) 


GRID(I) 


IGEOM 


TWOD 


LXTPLT 

LYTPLT 


Description 

Roberts  grid  transformation  sentinel  for  - 
direction.  , 

IGRID(l)  = 0.  Equally  spaced  grid  points. 

IGRID(i)  = 1.  One  boundary  layer  at  XGMIN(l) 
of  fractional  thickness  GRID(l). 

IGRID(l)  « 2.  Boundary  layers  at  both  XGyiN(l) 
and  XGMAX(I). 

IGRID(l)  = 3.  One  boundary  layer  at  XGyAX(l). 

Estimate  of  boundary  layer  thickness  for 
Roberts  grid  transformation  in  Xj  - direction; 
fraction  of  XGMAX(l)  - XGMIN(l). 

Geometry  option  sentinel 

IGEOM  = ± 1.  Carstesian  coordinate  system, 
IGEOM  = ± 2.  Cylindrical  polar  coordinate  sys- 
tem. A negative  value  of  IGECW  results  in 
printout  of  the  metric  scale  factors  and  their 
derivatives.  Default  value  is  1. 

Logical  variable  for  two-dimensional  mode  of 
operation  of  the  MINT  code. 

TWOD  = .FAIBE.  gives  the  3-D  mode  (defa'ult). 
TWOD  = .TRUE,  gives  the  2-D  mode. 

X^  grid  point  location  for  time  history 
printer  plots.  NXl  s LXTPLT  s NX2. 

X2  grid  point  location  for  time  history 
printer  plots.  NYl  s LYTPLT  s NY2. 


A L. 
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Namelist 

variable 

LZTPLT 

LXSPLT 

LYSPLT 

LZSPLT 

LXTPL2 

LYTPL2 

LZTPI2 

CLENG 

IDMPX 

IDMPY 

IDMPZ 

NXIN 

NYIN 

NZIN 


or 

name  , Description 

grid  point  location  for  time  history 
printer  plots.  NZl  s LZTPLT  « NZ2. 

grid  point  location  for  spatial  printer 
plots. 

Xg  grid  point  location  for  spatial  printer 
plots . 

X^  grid  point  location  for  spatial  printer 
plots . 

Not  used. 

Not  used. 

Not  used. 

Reference  length  (meters).  Also  appears  in 
$READ3  namelist. 

X2  grid  point  location  for  matrix  inversion  sub 
routine  error  check  printout  (MGAUSS  and  GAUSS) 

Xg  grid  point  location  for  matrix  inversion 
printout . 

X2  grid  point  location  for  matrix  inversion 
printout . 

X, -dimension  for  dynamic  storage  allocation  in 
the  dependent  variable  array  (AC).  Default  is 
NUMDX+2. 

X2-dlmension  for  dynamic  storage  allocation  in 
the  dependent  variable  array  (AC).  Default  is 
NUMDY+2. 

X- -dimens ion  for  the  dependent  variable  array. 
NZIN  = NUMDZ+2. 
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Description 


Total  number  of  grid  points  for  AC  array  in  main 
memory  ; set  by  program. 

Douglas-Gunn  subtraction  flag;  set  by  program. 
Not  used. 

Number  of  grid  points  in  difference  molecule; 
must  be  3. 

Not  used. 

Set  by  program  for  ADI  control. 

Not  used. 


Not  used. 


ICDC  = 0 for  sector  addressable  disk  usage  on 
UNIVAC . 

ICDC  = 1 for  record  addressable  disk  usage  on 
CDC. 

INCUT  = 0.  All  dependent  variables  stored  in 
main  memory.  No  mass  storage  required. 

INCUT  = 1.  Mass  storage  required  for  some  de- 
pendent variables . l/c  unit  number  NUNIT  = 11 
defined  in  main  program  DAL  is  used  for  this 
purpose. 

Parameters  set  by  program  for  control  of  mass 
storage  input/output  of  dependent  variables 

Reference  quantities 

Reynold's  number  based  on  reference  quantities; 
Re  = (calculated). 

Prandtl  number;  default  Pr  = C.73. 
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Namelist  or 

variable  name  Description 


SCMDT 


Schmidt  number;  default  Sc  = 0.73. 


PREFF 

SCEFF 

BOLTZ 

WREF 

DENSR 

HREF 

IREF 

PREF 

CPREF 

ZREF 

WCREF 

DREF 

QRREF 

VISCR 


Effective  t\irbulent  Prandtl  number;  default 

= 0.90. 


Effective  turbulent  Schmidt  nxjmber;  default 
®°eff  ^ 

Boltzmann  number;  Bo  = (calculated). 

Reference  velocity;  m/sec. 

3 

Reference  density;  default  = I.185  hg/m  . 

Reference  enthalpy;  h^^  = Cp  (calculated). 

Reference  temperature;  default  = 300°K. 

Reference  pressure;  p^  = 10  p^^R  T^/Wp  (calcu- 
lated) ® 


Reference  specific  heat;  defaxilt  c = 1004.8 
joule/kg- -K. 


Reference  reciprocal  molecular  weight;  Zjj  = 
lOvWjj  moles/kg. 


Reference  molecular  weight;  default  = 28.966 
gm/mole. 


Reference  binary  diffusion  coefficient;  not  used 
at  present. 

Reference  radiation  energy  flux;  joule/m^- 
sec. 

Reference  dynamic  viscosity;  default  = I.816 
X 10“5  kg/m- sec. 


i 
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GAMMA 

RGAS 

CLENG 

PNOM 

LAMFLO 

VISCL 

CTWO 


AVISCl, 
AVISC2, 
AVISO  3 


IH3TAG 


Ratio  of  specific  heats,  default  y = Used 

only  in  Mach  number  calculation  for  output. 

universal  gas  constant,  R '=  8.31^3  joule/°K- 
mole . ® 


Reference  length;  L,  meters. 

Nondimens ional  nominal  reference  pressure,  de- 
fault PNOM  =1.0. 

Laminar  or  turbulent  flow  sentinel. 

LAMFLO  = 0.  TurbxLLent  flow  (default). 

LAMFLO  = 1.  Laminar  flow. 

\ ' 

Nondimens ional  laminar  viscosity  at  standard 
conditions,  VISCL  = 1.0. 

Ratio  of  second  coefficient  of  viscosity  to  first 
coefficient  of  viscosity;  defa\ilt  is  2/3. 

Artificial  viscosity  coefficient  for  X -,  X -, 
X^-directions , respectively.  This  parameter 
equals  l/Re^,  where  Re^  is  the  computational- 
mesh  Reynolds  number  above  which  artificial 
viscosity  is  required  to  obtain  smooth  solutions. 
Default  value  for  each  is  0.5. 

IHSTAG  = 0.  Static  enthalpy  option.  Default. 

mSTAG  = 1.  Stagnation  enthalpy  option. 

IHSTAG  = 2.  Constant  stagnation  enthalpy  option. 
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1. 


Namelist  or 
variable  name 


Description 


IWFDIF 


Velocity  wall  function  differencing  flag. 
Not  used  when  IWVISC  = 3- 


IWFDIF  = 1.  Two-point  one-sided  differencing 
of  9u/dy. 

IWFDIF  = 2.  Central  differencing  of  Su/Sy. 
Default . 


IWVISC 


Turbulent  viscosity  wall  function  flag. 


IWVISC  = 1.  Wall  viscosity  consistent  with  wall 
slip  velocity. 

IWVISC  = 2.  Match  viscosity  derivative 
with  wall  function  value. 


YPL0, 

YPHI 


CLAMDA 


IWVISC  = 3.  Match  both  viscosity  derivative 
dUij,/Sy  and  velocity  derivative  Su/3y  using 
generalized  wall  function  formulation  (Ref.  3). 
Default . 

Lower  and  upper  limits  of  y"*"  for  which  van  Driest 
wall  function  is  employed.  Not  used  when 
IWVISC  = 3.  Default  YPLO  is  2.0,  YPHI  is  100.0. 

Constant  in  Lilley’s  mixing  length  model  (Ref.  3). 
Default  is  0.08. 


CWIIM 


Constant  in  Williamson's  mixing  length  model 
(Ref.  3).  Default  is  O.lU. 


FLILY 


Fraction  of  maximum  velocity  which  defines 
shear  layer  thickness  in  Lilley's  mixing  length 
model.  Default  is  0.0^. 


Von  Karmii  const. ant.  Default  is  O.k, 

Logarithmic  wall  function  constant,  (Ref.  3). 
Default  is  l.O/CVK. 


3.  Gibeling,  H.  J. , H.  McDonald,  and  W.  R.  Briley:  Development  of  a Three- 

Dimensional  Combustor  Flow  Analysis  . AFAPL-TR-75“59>  Volume  II,  October  1976. 
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Namelist  or 
variable  name 


Description 


I 

I 


1 


CLW2 

CTRBMU 

CTR31 

CTRB2 

SIGK 

SIGE 

FSTKE 

SIGMA 


Logarithmic  wall  function  constant,  C2  (Ref.  3). 
Default  is  5.1 


Two-equation  turbulence  model  constant,  C , 
(Ref.  3).  Default  is  O.09.  ^ 

Two-equation  turbulence  model  constant,  Ct , 
(Ref.  3).  Default  is  I.55. 

Two-equation  turbulence  model  constant,  Cg, 
(Ref.  3).  Defa\J.t  is  2.0. 


Two-equation  turbulence  model  constant,  a^, 
(Ref.  3).  Default  is  1.0. 


Two-equation  turbulence  model  constant,  c , 
(Ref.  3).  Default  is  1.3.  ® 


Free  stream  turbulence  kinetic  energy  level; 
fraction  of  0.5  Default  is  0.03. 

St  e fan -Bolt  znann  constant. 


RADLEN 

Characteristic  frequency-averaged  length  scale 
for  radiative  absorption. 

EWRAD 

Radiative  emissivity  of  solid  walls.  Default 
is  1.0  (black  body). 

$READ^ 

Miscellaneous  program  control  parameters 

NT 

Number  of  time  steps  to  be  taken.  Default  is 

NT  = 5. 


DT 


I 


Initial  nondimensional  time  step;  if  DT  < 0.0 
then  DT  is  set  to  |DT(  * CFL.  CFL  is  the  explicit 
stability  limit  time  step. 

Minimum  nondimensional  time  step  for  this  run. 

If  DTMIN  <0.0  then  DTMIN  is  set  to  [DIMIN]*  CFL. 


lU 


i 


DTMIN 


Namelist  or 
variable  name 

DTMAX 

IDTADJ 


NTSTEP 

NTCYST 

SSEPS 

DELDT 

PCNTl 

PCNT2 


Maximum  nondimensional  time  step  for  this  run. 

If  DTMAX<0.0,  then  DTMAX  is  set  to  | DTMAX  | <*CFL. 

Time  step  control  sentinel. 

IDTADJ  = 0.  Constant  DT  is  used  for  the  entire 
run. 

IDTADJ  = 1.  Time  step  is  adjusted  by  the  factor 
DELDT  based  on  specified  changes  PCNTl  and 
PCNT2 . Default . 

IDTADJ  = 2.  Time  step  is  cycled  between  the 
specified  values  DTMIN  and  DTMAX  using  an 
acceleration  parameter  concept.  A sequence  of 
NTSTEP  time  steps  is  used  under  this  option. 

Number  of  time  steps  for  cycling  option  used 
for  accelerating  convergence,  see  IDTADJ  des- 
cription. Default  is  NTSTEP  = 3. 

Time  step  number  during  current  run  at  which 
time  step  cycling  option  (IDTADJ  =2)  is  to  be 
initiated. 

Steady-state  convergence  criteria.  Default  is 
SSEPS  = 0.001. 

Factor  to  increase  or  decrease  time  step  under 
IDTADJ  = 1 option.  Default  is  DELDT  = 1.25. 

If  maximum  relative  change  in  any  flow  variable 
is  less  than  PCNTl,  the  time  step  is  increased 
to  DELDT  * DT  (IDTADJ  = 1 option).  Default  is 
PCNTl  = O.OU. 

If  maximum  relative  change  in  any  flow  variable 
is  greater  than  PCNT2,  the  time  step  is  de- 
creased to  DT/DEIDT  (IDTADJ  = 1 option).  De- 
fault is  PCNT2  = 0.06. 
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Namelist  or 
variable  name 


Description 


TTINE 


ESTIME 


IPRINT 

ITEST 

IPLOT 

ITPLT 

NTREST 

IDUMPl 


Nondimens ional  value  of  problem  time  at  start 
of  run.  For  steady  computations  this  value 
enters  only  into  the  pseudo-kinetic  hydrocarbon 
chemistry  rate  constant. 

Estimated  computer  run  time  in  minutes.  If  the 
run  time  measured  from  the  start  of  execution  of 
this  program  exceeds  ESTIME  at  the  end  of  any 
tine  step,  the  program  will  terminate  execution 
after  a normal  final  print.  Default  is  ESTIME 
= 600.0  minutes. 

Complete  flow  field  printouts  are  provided  every 
IPRINT  time  steps.  Default  is  IPRINT  = 100 

Steady  state  test  is  performed  every  ITEST  tine 
steps.  Default  is  ITEST  = 1. 

Not  used  at  present. 

Number  of  time  steps  for  which  selected  time 
plot  data  are  saved  on  restart.  Set  by  program. 

Time  step  interval  at  which  restart  dirr.ps  are 
taken.  See  IREPUN  description  under  $F.EAD1 
namelist. 

Supplementary  print  flag. 

IDUMPl  = 1.  Print  initial  flow  field. 

IDUMPl  =2.  No  initial  printout. 

IDUMPl  =-l.  Print  initial  flow  field  and 
solution  variables  at  the  end  of  each  ADI 
step  for  the  first  time  step. 

IDUMPl  = 0.  Print  solution  variables  at  end 
of  each  ADI  step  for  first  time  step. 

Default  is  IDUMPl  = 1. 
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Namelist  or 
variable  name 


Description 


MGPRJJT 

NAMLST 

TEXIT 

TSTOP 

MEQS 

MCPLD 

NEQN(IP,IADI) 

NEQN(IP+1,IADI) 

NPADI 


Matrix  inversion  routine  error-check  printout 
is  provided  every  MGPRNT  time  steps.  This  print 
is  automatically  provided  at  the  first  time  step. 
Default  is  MGPRNT  = NT. 

Sentinel  for  optional  namelist  input  in  subrou- 
tine BVIV.  NAMLST  =0:  No  additional  input 

read. 

If  the  computer  run  time  remaining  at  the  end  of 
any  time  step  is  less  than  TEXIT,  the  program 
will  terminate  execution  after  a normal  final 
print.  Default  is  TEXIT  = 20.0  seconds. 

If  the  nondimensional  problem  time  TTIME  exceeds 
TSTOP,  the  program  will  terminate  execution  in 
the  normal  manner.  Default  is  TSTOP  = 10^*^. 

The  total  number  of  equations  to  be  solved  during 
the  flow  field  calculation  procedure  (NOX  chemis- 
try is  not  included).  The  actual  equations  are 
specified  via  the  HftNUM  array.. 

The  total  number  of  equations  to  be  solved 
coupled  during  each  ADI  sweep.  Used  only  when 
NPADI  = 1. 

The  number  of  coupled*  equations  to  be  solved 
during  sequence  number  IP  of  sweep  lADI. 

The  number  of  decoupled*  equations  to  be  solved 
during  sequence  nximber  IP  of  sweep  lADI. 

It 

The  number  of  sequences  of  coupled  and 
decoupled  equations  to  be  solved  during  each 
ADI  step. 


*The  notation  "coupled"  and  "decoupled"  is  used  here  to  indicate  coupling 
between  the  finite  difference  equations  during  any  ADI  step. 


Namelist  or 
variable  name 


Description 


INPEQS 

rEe,NU!/(j) 

IVARNO(I) 

ISOLVE(I) 


ICPLD(I,IADI) 


* 'lee  previous 


Not  used. 

Array  containing  the  external  equation  numbers 
(see  Table  JI ) of  the  equations  to  be  solved; 

1 S j s MB5.S. 

The  index  of  the  dependent  variable  associated 
with  external  equation  I = lEQMt'(j).  See 
Table  II . 

Flag  which  specifies  if  the  external  equation 
I = lEQjrjM(j)  is  to  be  solved  during  the  present 
run,  or  at  a later  time  via  a restart. 

ISOLVE(I)  =0.  Do  not  solve  equation. 

ISOLVE(i)  = 1.  Solve  equation. 

ISOLVE(I)  s 2.  Solve  equation  at  a later  time, 
by  specifying  ISOLVE  = 1 in  the  restart  mode. 
Devault  values  are  given  in  Table  III  . 

Coupled  equation  sentinel  for  external  equation 
n'umber  I on  ADI  sweep  lADI  (lADI  = 1,2,3).  De- 
fault values  are  set  in  main  program  DAL.  See 
Table  HI.  The  meaning  of  these  values  is  as 
follows,  where  IP  = ICPLD (l,IADl) . 

I s I?  « 10.  Equation  number  I is  coupled* 
during  sequence  number  IP  of  sweep  lADI. 

II  S IP  s 20.  Equation  number  I is  decoupled* 
d'uring  sequence  number  (IP  - 10)  of  sweep  lADI. 

21  S IP  S 30.  Equation  number  I is  independent 
of  flow  field  during  sequence  number  (I?  - 20) 
of  sweep  lADI. 

IP  = 0.  Equation  number  I is  not  solved  during 
sweep  lADI. 

footnote. 
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Namelist  or 
variable  name 


JBQN(I,IADI) 


JVAR(I,IADI) 


Description 

External  equation  number  of  the  Ith  equation 
which  is  solved  during  sweep  lADI.  Set  by 
program. 

Variable  number  associated  with  the  Ith 
equation  which  is  solved  during  sweep  lADI. 
Set  by  program. 


NVDMP  The  number  of  dependent  variables  in  the  AC 

array  to  be  written  on  the  restart  d\jmp.  Only 
a single  time  level  is  written  for  these 
variables.  Default  is  NVDMP  = MB9.S. 


NVDMPl,  Indices  for  dumping  the  remainder  of  AC  array 

NVDMP2  req\iired  for  a restart.  Set  by  program. 

NVDMPl  = NVDMP  + 1;  NVDMP2  = last  variable 
index  being  used  by  program. 

The  following  variables  are  integer  indices  used  in  referencing  the  correct 
component  of  the  AC  array  in  the  program.  These  variables  are  stored  at  two 
time  levels . 


U 

V 

w 

R 

H 

MF 

MFP 

TKE 

TED 


u - velocity  component  index 
V - velocity  component  index 
w - velocity  component  index 
Density  index 

Stagnation  or  static  enthalpy  index 

Uhburned  fuel  mass  fraction  (m^)  index 

Mixture  fraction  ({)  index 

Turbulence  kinetic  energy  index 

Turbulence  kinetic  energy  dissipation  rate 
index 
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Namelist  or 
Variable  name 

VS 

NO 

NPARTl 

NPART2 

The  following 
stored  in  the 

T,  TL 
ME, MEL 

CH,CHL 

ZN,ZNL 

D,DL 

P,PL 

VAP,VAPL 

CK1,CK1L 

CK2,CK2L 

CKMjCKMJ. 

CKH,CKHL 

CKP,CKPL 

NOljNOlL 

N02,N02L 


Description 

Effective  viscosity  index 

Nitrogen  oxide  mass  fraction  index 
Lower  particle  fraction  equation  nijmber 

Upper  particle  fraction  equation  number 

variables  are  integer  indices  used  in  referencing  q'oantitie 
AC  array  at  only  a single  time  level. 

Temperature  index 

Equilibrium  hydrocarbon  fuel  mass  fraction 
(m^g)  index 

Reciprocal  specific  heat  index  (c„  = l/c  ) 

Reciprocal  molecular  weight  (Z)  index 

Velocity  devergence  index 

Static  pressure  difference  (p-p  ) index 

nom 

Droplet  vaporization  source  term  index 
I^drocarbon  chemistry  function  index 
}]^drocarbon  chemistry  function  index 
Hydrocarbon  chemistry  function  index 
Hydrocarbon  chemistry  function  index 
Hydrocarbon  chemistry  function  index 
Nitrogen  oxide  chemistry  function  index 
Nitrogen  oxide  chemistry  function  index 
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Namelist  or 
Variable  name 


Description 


N03,N03L 

NOU,NOl+L 

DLIMIT 

IDIFF(IADI 

NCOMP 

NPART 

PSTATl 

RHOl 

PDROP 


Nitrogen  oxide  chemistry  function  index 
Nitrogen  oxide  chemistry  function  index 
Not  used 

lER)  Sentinel  for  optional  nonlinear  differencing 

formulas. 

IDIFF  = 0,  Employ  conventional  central  differ- 
ence formulas.  Default. 

IDIFF  = 1.  Employ  nonlinear  difference  formulas 
for  ADI  sweep  number  lADI  in  equation  number  lEQ. 

Number  of  components  for  multicomponent  flow. 
NCOMP  = 1.  Single  component  flow  of  air. 

NCOMP  = 2.  Two-component  flow  of  methane  and 
air  (default).  Other  species  may  be  treated  by 
reading  in  the  appropriate  thermodynamic  data  in 
JANNAF  format  into  the  CG  array  in  the  $CHEMD 
namelist. 

Number  of  particle  fraction  classes  for  droplet 
vaporization  model.  NPART  s 4.  Default  is 
NPART  = 0,  which  means  that  droplet  vaporization 
will  not  be  considered.  Note  that  program  di- 
mension NVAR  in  the  AC  array  must  be  increased 
by  NPART  for  storage  of  droplet  particle  frac- 
tions . 

Average  nondimens ional  inlet  static  pressure. 
Default  is  PSTATl  =1.0. 

Not  used. 

Pressure  drop  parameter.  The  specified  exit 
static  pressure  is  set  to  PSTAT  = PDROP  * PSTATl. 
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Namelist  or 
variable  name 


Description 


R>!FLW1 

CMFLW(J) 

IPGRAD 


LZPRNT(LZ) 


IVARPR(IV) 


Not  used. 

Not  used. 

Pressure  differencing  parameter. 

IPGRAD  = 0.  Central  differencing  of  pressure 
in  all  coordinate  directions.  Default. 

IPGRAD  = 1.  C>ne-sided  differencing  in  X-, - 
and  X2  directions. 

IPGRAD  = 2.  Nonlinear  differencing  in  X-,  and 
Xg  directions. 

Optional  print  control  flag. 

LZPRNT  = 0.  Suppress  all  flow  field  printout 
at  axial  station  number  LZ. 

LZPRNT  = 1.  Normal  flow  field  printout  at  axial 
station  LZ.  Default  is  1. 

Optional  print  control  flag  for  dependent  vari- 
able n’umber  IV.  (Table  II  ). 

IVARPR  = 0.  Suppress  printout  of  variable  IV. 

IVARPR  = 1.  Normal  printout  of  variable  TV. 
Defaijlt. 


$READ3 


Boundary  conditions  and  wall-port  specification 


IBOUND(rBC) 


Boundary-type  indicator  for  surface  IBC  (See 

Figs.  3, 

u). 

See  Table  V for  default  values . 

IBOUND  = 

1. 

Inlet  or  exit  plane  (no  walls). 

IBOUND  = 

2. 

Plane  of  symmetry. 

IBOUND  = 

3. 

Free  surface. 

IBOUND  = 

1+. 

.Wall -function  surface*. 

IBOUND  = 

5. 

Wall-function  surface*  with  ports 

IBOUND  = 

6. 

No  slip  wall. 

*A  wall-function  surface  is  a wall  at  which  a slip  velocity  is  determined 
using  the  law-of-the-wall. 


t 
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Namelist  or 
variable  name 


Description 


naBC(lBC,l)  Boundary  condition  specification  for  external 

eqiiation  number  I at  svirface  IBC  (See  Figs.  3,14). 
See  Table  IV  for  default  values. 

Standard  boundary  conditions  (no  testing  for 
ports ) : 

lEftBC  = 1.  Function  boundary  condition  (depen- 
dent variable  only) . 

lEftBC  = 2.  First  derivative  of  dependent  vari- 
able eqiials  0.0. 

HftBC  = 3-  Second  derivative  of  dependent  vari- 
able equals  0.0. 

Special  boundary  conditions  (test  for  ports  if 
IBOUND(IBC)  = 5): 

Continuity  equation  (1=4) 

HftBC  =4.  At  a wall,  apply  time-dependent  con- 
tinuity equation.  At  a port,  specify  density  or 
apply  a pressure  bovmdary  condition,  as  deter- 
mined by  ICONBC(IBC),  see  below. 

Energy  equation  (I  = ^) 

HftBC  s 4.  At  a wall,  apply  adiabatic  wall 
condition. 

HfiBC  =4.  At  the  inlet,  specify  static  tem- 
perature (see  initial  flow  field  input  descrip- 
tion). Cannot  be  used  when  density  is  speci- 
fied at  the  inlet. 

IBiBC  = Not  used. 

URBC  = 6.  dp/Sn  = 0 where  p is  static 
pressure  and  n is  normal  to  surface  IBC. 

UftBC  = 7.  d^p/dn^  = 0 (recommended). 

Momentum  equations  (l  = 1,2,3)  wall  boundary 
conditions; 

HftBC  = 4.  Wall  function  boundary  condition  (no 
ports  permitted). 

HftBC  4 5.  Wall  function  boundary  condition  at 
wall  grid  points. 
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Namelist  or 
variable  name 


Description 


IEQBC(IBC,I) 

continued 


yoment'-E  equations  (l  = 1,2,3)  and  species  equa- 
tions (I  = 6 to  12)  port  boundary  conditions. 
IES3C  = 5.  Function  boundary  condition,  i.e., 
specified  velocity,  mass  fraction  or  particle 
fraction  profiles . 

IBQ3C  = 6.  S^j/Bn  = 0 
IE0,BC  = 7.  = 0 

Species  equations  (I  = 6 to  12)  wall  boundary 
condition 


HftBC  a 4.  First  derivative  of  mass  fraction  or 
particle  fraction  equals  zero. 


Radiation  transport  equations  (l  = 13  to  1^) 

IB93C  =4.  At  a wall,  the  specified  wall  temp- 
erature (T„)  boundary  condition  is  applied.  At 
present  = 300°K  is  employed  and  the  wall 
emissivity  is  taken  as  e =1.  At  a port,  out- 
going radiation  escapes  and  the  appropriate 
boundary  condition  requires  a specified  in- 
coming radiative  flux,  q^^  At  present  q.^  = 
0.0  is  employed. 


Turbulence  kinetic  energy  equation  (l  = 16) 
lEQBC  =4.  At  a wall,  linearly  extrapolate 
k from  interior  points  using  h^'k/'dr?  = 0. 
lEQBC  = 5-  Inflow  port  designation,  specified  k. 
lEQBC  = 7.  Outflow  port  designation,  d^k/'dn^  = 0. 


Turbulence  kinetic  ener^  dissipation  rate 
equation  (I  = 17) 

IBOBC  24.  At  a wall,  c is  specified  consistent 
with  wall  function  formulation  at  one  grid  point 
away  from  wall. 

lES^C  = 5.  Inflow  port  designation,  specified  s. 
IEQ3C  = 7.  Outflow  port  designation,  ?^e/dn^  = 0. 
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Namelist  or 
variable  name 


Description 


I 


I 


i 


i 


ICONBC(IBC) 


NPORTS 

LCPORT(NP) 

n>0Rl(lUJ,IG2,NP) 


IP0R2(ILU,IG1,NP) 


Continuity  equation  port  boundary  condition  in- 
dicator for  surface  IBC,  See  Table  IV. 


ICONBC  = 1.  Specified  density  at  port  (re- 
commended at  inlet  or  injection  ports). 

ICONBC  = 2.  Specified  static  pressure  at 
port  (recommended  at  exit  port) 

ICONBC  = 3.  Specified  stagnation  pressure  at 
port. 


ICONBC  s U, 
ICONBC  = 5. 
ICONBC  = 6. 


Specified  Sp/Sn  at  port. 
3p/3n  = 0.0  at  port. 
b p/dn^  = 0.0  at  port. 


Total  number  of  ports  on  all  boundaries . 
NPORTS  S 6. 


Must  be  set  equal  to  the  surface  niimber  IBC 
(Figs.  3 and  U)  containing  port  number  NP. 


Port  coordinate  specification.  When  HU 
equals  1 (or  2)  IPORl  designates  the  lower  (or 
upper)  grid  point  number  (for  port  number  NP) 
in  the  Xi  - coordinate  direction  at  grid  point 
number  IG2  in  the  X2  - coordinate  direction. 

See  footnote  below. 

Port  coordinate  specification.  When  ILU 
equals  1 (or  2)  IP0R2  designates  the  lower 
(or  upper)  gjid  point  number  (for  port  number 
NP)  in  the  X^  - coordinate  direction  at  grid 
point  number  IGl  in  the  X^  - coordinate  direc- 
tion. See  footnote  below. + The  specification 
of  IP0R2  must  be  consistent  with  that  of  IPORl, 
and  it  is  req\iired  to  define  a general  port  ge- 
ometry euid  to  avoid  excessive  internal  testing 
in  the  program. 


Note  on  port  coordinate  specification.  The  actual  coordinate  directions 
corresponding  to  X^  and  X^  are  determined  by  the  surface  (IBC)  which  contains 
the  port  NP,  as  Indicated  in  the  following  table: 
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* 

* 

IBC 

"2 

1,2 

- 

1 

1 

5,6 

^1 

Also,  the  grid  point  numbers  specified  in  the  IPORl  and  IP0R2  arrays  are 
assumed  to  coincide  with  the  wall  grid  point  closest  to  the  port  under  con- 
sideration. If  the  computational  boundary  coincides  with,  a line  of  symmetry 
through  the  port,  then  the  values  of  IPORl  (or  IP0R2)  defining  that  pert 
boundary  must  be  set  to  zero. 


$CHEMD 


Chemistry  and  droplet  data 


ICHEMR 


Hyrdocarbon  chemistry  sentinel 
ICHEMR  = -1.  No  chemical  reactions, 
variable  specific  heat  (c  ) using  JANNAF 
data.  ^ 


ICHEMR  = 0.  No  chemical  reactions,  con- 
stant specific  heat  (c^).  See  CCP  in- 
put description  below.* 

ICHEMR  = 1.  Pseudo-kinetic  hydrocarbon 
chemistry,  variable  specific  heats  (c  ) 
using  JANNAF  data. 


JCHEM  Chemistry  flag  set  in  subroutine  TEMPN 

ATC,ATH  Hydrocarbon  fuel  specification,  C H , 

where  n = ATC,  m = ATH.  Defaults'^are 
ATC  =1.0,  ATH  = h.O  (methane). 

OXMF  Mass  fraction  of  ox-/gen  (Og)  in  the  ox- 

idizer (normally  air).  OXMF  = 0.2322 
for  air. 

TFLOOR  Nondimensional  minimum  temperature  per- 

mitted with  chemistry.  Default  is 
TFLOOR  = 300°K/TREF. 
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Namelist  or 
variatole  name 


CHFM(J) 

CMW(J) 

CMWINV(J) 

CCP(J) 


I 


I 


IFUEL 


EPEC 


TCHEMl, 

TCHIW2, 

TCHEM3 


I TSCAL1,TSCAI2 


Description 

Heat  of  formation  of  component  J (kcni/ 
mole) 

J = l:  methane  (CH4),  CHFM  = 17.895 
kcal/mole . 

J = 2:  air,  CHFM  = 0.0. 

Molecular  weight  of  component  J (gm/mole). 
J = 1:  methane,  CMW  = l6.043. 

J = 2;  air,  CMW  = 28.966. 

Inverse  of  CMW(j),  calculated. 

Specific  heat  of  component  J (Joule/ 
kg-°K)  for  ICHEMR  = 0 option. 

J = 1:  methane,  CCP  = l88t.O  Joule/ 

kg-°K. 

J = 2:  air,  CCP  = 1004.8  Joule/kg-°K. 

% 

Fuel  species  index  for  JANNAF  tabular 
data.  Default  is  IFUEL  = 12,  see  BLOCK 
DATA  subprogram  (CHIMDATA)  and  CG  des- 
cription below. 


Total  number  of  species  in  the  3ANNAF 
table  (CG  array,  see  below).  EPBC  = 15. 

TCHEML  is  the  nondimensional  start 
time  for  pseudo-kinetic  hydrocarbon 
chemistry.  If  TCHEMl  s TTIME  s TCHEM2, 
use  TSCALl  as  the  iC  chemistry  kinetic 
time  scale.  If  TCHE^C  < TTIME  s TCHEM3, 
decrease  the  kinetic  time  scale  from 
TSCALl  to  TSCAI2  vising  a cosine  function 
of  time.  If  TTIME  > TCHEM3,  use  TSCAI2 
as  the  1C  chemistry  kinetic  time  scale. 

Nondimensional  time  scales  for  pseudo- 
kinetic  HC  chemistry  calculation. 
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Namelist  or 
variable  name 


IRC  ON 


I 

1 


ZMNRAT 

ZMXRAT 

LZRATl 

LZRAT2 

FRATE 

EPST 

ITERT 


Description 

Rate  constant  (or  time  scale)  control 
sentinel  for  pseado-kinetic  KC  chemistry/' 
calc’alation.  Determines  the  tjrpe  of  ig- 
nition delay  used  in  the  HC  chemistry 
calculation. 

IRCON  = 0.  No  delay;  i.e.,  rate  constant 
independent  of  distance  from  inlet  rort. 

IRCON  = 1.  Linear  variation  of  rate 
constant  r^  as  a function  of  distance 
from  inlet  port,  where  (TSCAL1}““^ 

Pc  s (TSCAI2)-!  for  ZMNRAT  ^ X3  s zMXRAT. 

IRCON  = 2.  Quadratic  variation  of  rate 
constant  as  a function  of  distance  from 
inlet  port. 

Nondimens ional  value  of  X^-coordinate 
such  that  the  rate  constant  equals 
(TSCALl)"^  in  the  range  ZMIN  ^ Xq  s 
ZMNRAT.  Default  is  ZMNRAT  = Z(3J. 

Nondimensional  value  of  X^-coordinate  at 
which  the  rate  constant  reaches  the  maxi- 
mum value  (TSCAI^)"-.  See  IRCON  discus- 
sion above.  Default  is  ZMXRAT  = ZMAX. 

Optional.  If  specified  then  ZMNRAT  is 
set  equal  to  Z( LZRATl). 

Optional.  If  specified  then  ZMXRAT  is 
set  equal  to  Z(LZRAT2). 

Rate  constant  mviltiplier.  Default  is 
FRATE  =1.0. 

Not  used  at  present. 

Not  used  at  present. 
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Namelist  or 
variable  name 


Description 

Nitrogen  oxide  chemistry  sentinel 


NOX 


NCK  =0.  Do  not  solve  NO  chemistry 
equat ion . 

NOX  = 1.  Solve  NO  chemistry  equation 
after  convergence  of  flow  field  solu- 
tion. This  option  is  not  operational. 

NOX  = 2.  Solve  NO  chemistry  equation 
using  flow  field  obtained  from  a 
restart  dump. 

NTNOX  Number  of  time  steps  for  NO  chemistry 

calculation.  Overrides  input  value 
of  NT. 

NOXDN  Output  unit  number  for  partial  equilibrium 

mole  fraction  data  calculated  during  the 
nitric  oxide  chemistry  initialization 
procedure  in  subroutine  CNOXI.  The 
partial  equilibrium  calculation  is 
fairly  time  consuming,  so  it  is  de- 
sirable to  save  the  necessary  mole 
fractions  on  a disk  or  tape  file.  De- 
fault is  NOXUN  = 21. 

NOREST  Nitric  oxide  chemistry  restart  option. 

See  NOXUN  above. 

NOREST  =0.  Do  not  write  partial 
equilibrium  mole  fraction  data  on  unit 
NOXUN. 

NOREST  = 1.  Write  partial  equilibrium 
mole  fraction  data  on  unit  NOXUN 
(default). 

NOREST  = 2.  Read  partial  equilibrium 
mole  fraction  data  from  unit  NOXUN. 
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Namelist  or 
variable  name 


FNOXST 


ANOX(J), 

BNOX(J), 

PNOX(J), 


Description 


The  initial  guess  for  nitric  oxide  lass 
fraction  is  set  equal  to  FNOXST  times 
the  value  obtained  from  the  partial 
equilibrium  calculation.  Defa’jlt  is 
FNOXST  = 0.01. 

Chemical  rate  constants  for  extended 
Zeldovich  nitric  oxide  reaction  mecha- 
nism. An  Arrenhius  form  is  assumed; 
kj  = A.T^J  exp  (-B.,/RT).  See  Table  V 
for  values  and  ujii^s. 


The  following  input  parameters  pertain  to  the  partial  equilibrium  calculation 
(subrouting  EftUIL)  needed  for  the  nitric  oxide  chemistry  initialization  pro- 
cedure (subroutine  CNOXi). 

FA.MAX  Maximum  fuel-to-air  mass  ratio  for  which 

the  equilibrium  calculation  will  be  per- 
formed. Default  is  FA.MAX  = l.OxlO"". 

EPSEJ^M  Convergence  criterion  for  equilibrium 

Newton-Paphson  iteration  procedure  (sub- 
routine EQUIL).  Default  is  EPSEQ>!  = 1.0 
X 10"  5 

EPSEXiT  Convergence  criterion  for  initial 

temperature  guess  (needed  by  EQUIL) 
iteration  procedure  (subroutine  SETl). 
Default  is  EPSEQT  = l.OxlO"^. 
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Namelist  or 
vaj-jable  name 

TMIN 

TMAX 

ICSiMX 

ntmaj: 

The  following 
RDROP(J) 

RHODR 

HDVAP 

HDCOMB 

BSTC 


Description 


Minimum  temperature  allowed  during 
EftUIL  iteration.  Default  is  TMIN  = 
290°K. 

Maxim'im  temperature  allowed  during  B^UIL 
iteration.  Default  is  TMAX  = 6000°K. 

Maximum  niimber  of  iterations  for  05,UIL 
iteration  procedure.  Default  is  ICE^l^K 

= 25. 

Maximum  number  of  iterations  for  3ET1 
iteration  procedure.  Default  is  NTMAX 
= 10. 


input  parameters  pertain  to  the  droplet  vaporization  model. 


Discrete  droplet  radii  (meters)  for  the 
particle  fraction  classes  (see  Fig.  5); 
1 s J S NPART  + 1 and  RDROP(K)  > 
RDR0P(K+1)  is  required.  Defa\iLt  wlues 
are  1.25x10-^,  l.OxlO"^,  0.75x10"^, 
0.5x10"^,  0.25xl0"^m,  respectively. 


Droplet  (liquid)  density,  a constant. 
Default  is  RHODR  = 640.7  kg/m^  (liquid 
CH^). 

Heat  of  vaporization  of  the  liquid  drop- 
lets. Default  is  HDVAP  = 0.31966x10^ 
Joule/kg  (liquid  CH4). 

Heat  of  reaction  for  stoichiometric  com- 
bustion of  the  fuel  with  O.^.  Default 
is  HDCOMB  = 50.2x10°  joule/kg  (CH^). 

Ratio  of  fuel  mass  to  oxygen  mass  for 
the  stoichiometric  combustion  reaction. 
Default  is  RSTC  = 0.25  (CHj^  + 2O2  - 
COg  + BHgO). 
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Namelist  or 
variable  name 


Description 


RHDVAP 


Inverse  of  HDVAP.  Calculated. 


CRDRI(J), 

CRDRP(J) 


Calculated  f’unctions  of  the  radii  RDROP: 
1 5 J s NPART  + 1. 


The  following  input  arrays  contain  the  chemical  species  thermodynamic  data. 
Polynomial  fits  to  the  JANNAF  data  have  been  made  over  cwo  temperatiire  ranges 
(300  - 1400°K  and  14oO-6000°K) . Seven  polynomial  coefficients  are  stored  for 
each  required  species  (see  Table  VI)  for  each  tempera t’ure  range. 


C G ( I , IS ) Polynomi al  coefficients  for  species 

number  IS . Temperature  range  1^00  - 
6000°K,  1 « I 5 7;  temperature  range 
300  - 1400°K,  8 « I 14.  The  format 
of  CG  is  implied  by  the  following  non- 
dimensional  polynomial  relationships, 
with  Cj  designating  CG( I, IS)  or  CG( 1*7, 
IS)  for  the  two  temperature  ranges  con- 
sidered; 

C^  = C.  + C,T  + C,T^+  CkT  + CcT 
^ 1 3 5 

R 

H = C^  + C2T  + C^^T'^  * C^T^  + C^T^  * C^ 
RT  2 3 4 5 T 


where  R is  the  universal  gas  constant, 
and  Cp  and  H are  the  molar  specific  heat 
and  molar  enthalpy,  respectively.  The 
temperature  T is  in  degrees  Kelvin,  which 
implies  the  units  for  the  coefficients 
Cj  (l  = 1 to  6).  The  coefficient  C is 
not  required  in  the  present  analysis. 

SMW(IS)  Molecular  weight  of  species  IS,  gm/mole. 

See  Table  VI. 

ACP(I,J)  Nondimen si on al  polynomial  coefficients 

for  fuel  (J  = 1)  and  air  (j  = 2)  calcu- 
lated from  CG  data. 


Namelist  or 
variable  name 

SHFM(IS) 


BETA(IS,JE) 


Description 

Nondimensional  heat  of  formation  of 
species  IS  calculated  from  CG  data. 

Atomic  stoichiometric  coefficients; 
the  number  of  gm-atoms  of  element 
JE  per  gm-mole  of  species  IS.  See 
Table  VII. 


Namelist  Input  For  Initial  Flow  Field 


This  input  is  read  in  subroutine  FLWFID  only  if  IP.SST  = 0,  and  following 
the  previous  namelist  $READ1  through  Namelist  $30'J!O  is  read  once, 

and  then  namelist  $BCDATA  is  read  at  each  - coordinate  location  (see  Fig. 
h)  for  which  LZRE.AD(LZ)  is  greater  than  zero,  1 S LZ  s NZ2,  and  for  either  a 
user  specified  number  of  X2  - coordinate  locations  (see  NTfREAD)  or  for  a 
specified  number  of  boundary  locations  (see  NBCRD). 

Namelist  or 

Variable  name  Description 

$B0UND  Initisil  Flow  Field  Setup 


LZREAD(LZ)  Input  sentinel  for  X,  - grid  point  number  LZ. 

LZREAD  =0.  No  input  will  be  read  for  this  LZ 
station;  program  sets  up  a default  flow  field 
based  on  inlet  profile  data. 

LZREAD  = 1.  Read  boundary  data  into  the  array 
AIE  for  this  LZ  station  (see  below'. 

LZREAD  = 2.  Read  entire  starting  solution  into 
the  array  AIE  for  this  LZ  station  (see  below). 
LZREAD  = 3*  Same  as  LZREAD  = 2,  except  the 
X^  - velocity  component  will  be  a power  law  pro- 
file. Intended  for  pipe  and  channel  flews. 


NYREAD(LZ) 


NBCRD(LZ) 


Number  of  X2  - coordinate  locations  for  this  LZ 
station  at  which  namelist  $BCDATA  will  be  read. 
Defa-ult  is  NEMDY  + 2.  Used  only  when  LZREAD  (LZ) 
i 2.  See  INY2  description  below. 

Number  of  boundary  locations  IBC  (1  to  h)  for  thi 
LZ  station  at  which  boundary  data  only  will  be 
read  via  namelist  $BCDATA.  Default  is  U.  Used 
only  when  LZREAD (LZ)  = 1.  See  INIBC  description 
below. 


IBPOWR 


Not  used. 


FNTURB,  RADIUS,  Parameters  used  for  defining  optional  power  law 

WINCL,  XJCL,  yjCL  velocity  profile  (LZREAD  =3).  FNTURB  is  the 

power  law  exponent  for  the  profile  u^/u  = 

(1  . where  n = FNTURB,  R = RADIuS, 

Upj^  = WHK;L,  and  r is  the  radieil  coordinate 
measured  from  the  location  (x^,:^)  = (XJCL,  YJCL) 


Namelist  or 
variable  name 


Description 


IVS0PT 


INVS 


VISCIN 

Nvmp 

lEQINP(lVAR) 


I 

ISPRED 


DECAY 


NJETS 


Turbulent  viscosity  option  used  for  computing 
initial  for  two-equation  turb\J.ence  model 
initialization. 

IVS0PT  = 2.  Code  utilizes  procedure  described  in 
Ref.  3. 

IVS0PT  = 3-  Turbulent  viscosities  are  supplied 
on  input  (see  INVS  below),  thereby  bypassing  the 
programmed  mixing  length  model  formulation. 

Default  is  2. 

Turbulent  viscosity  may  be  loaded  into  AIE(LG,INVS) 
in  namelist  $BCDATA  when  IVS0PT  = 3 option  is 
being  used.  Default  INVS  = 7. 

Nondimen sional  dynamic  viscosity  for  laminar 
flow.  Default  is  VISCIN  = 1.0 

Number  of  variables  to  be  read  into  the  AIE  array 
described  below.  Default  NVINP  = 7. 

Array  to  relate  the  input  variable  number  IVAR 
(1  < IVAR  < NVINP)  to  the  external  equation  num- 
ber lEft  (Table  II).  The  default  values  are  given 
in  Table  VIII. 

ISPRED  =0.  Do  not  spread  and  decay  inlet 
velocity  profile. 

ISPRED  = 1.  Spread  and  decay  inlet  velocity 
profile  based  on  free  jet  similarity  solution. 
Default  is  ISPRED  = 1 if  there  is  an  inlet  port. 

Velocity  decay  exponent  for  two-dimensional 
jets  (DECAY  =0.5)  and  circular  jets  (DECAY  = 

1.0). 

Number  of  inlet  port  jets.  Default  NJETS  = 1. 


35 


Namelist  or 
variable  name 

LXPC(IJ), 

LYPC(IJ) 


XJZERO(IJ) 

yjZERO(lJ) 

LXSPRD(I,IJ) 

LySPRD(l,IJ) 

LZSPRD(I,IJ) 


$BCDATA 

AIE(LG,IVAR) 


INYl 


INY2 


INIBC 


NVINP, 

miNP(lVAR) 


Description 


The  X]_  and  X2  grid  point  coordinates,  respec- 
tively, for  the  port  center  of  inlet  port  number 
IJ.  Must  be  specified  unless  the  actual 
physical  coordinates  XJZ3P0(lJ)  and  yjZSHO'lJ) 
are  specified  instead. 

The  and  X^  physical  coordinates,  respectively, 
for  the  port  center  of  inlet  port  number  IJ. 

The  X]_,  Xp,  and  X^  grid  point  coordinates, 
respectively,  that  determine  the  points  in  the 
(Xj_,  X-)  plane  and  in  the  (X,,  X^)  plane  to 
which  the  inlet  jet  will  be  spread.  The  indices 
1=1  and  1=2  designate  the  lower  and  upper 
coordinates  of  the  jet  emerging  from  port 
nvnnber  IJ. 

Flow  field  data  input 

Nondimensional  value  of  variable  number  IVAR 
(Table  VIII)  at  grid  point  (LX,LY)  for  the  current 
LZ  station,  where  LG  = IX  and  LY  is  determined 
by  INYl  and  INY2  as  described  below.  When 
boundary  data  is  being  read  (LZREAD  =1),  LG  = LY 
if  INIBC  = 1 or  2 and  LG  s LX  if  inibC  = 3 or  U. 

First  LY  (X2  coordinate)  grid  point  number  for 
the  data  being  loaded  into  the  AIE  array  during 
the  present  $BCDATA  namelist  read.  Set  by  code. 

Last  LY  (X2  coordinate)  grid  point  number  for  the 
data  being  loaded  into  the  AIE  array  during  the 
present  $BCDATA  namelist  read.  Must  be  speci- 
fied. 

IBC  number  (Fig.  3 ) for  the  boundary  data  being 
loaded  into  the  AIE  array  during  the  present 
$BCDATA  namelist  read.  Used  only  when 
LZREAD(LZ)  = 1. 

If  specified  these  values  will  override  those  set 
in  namelist  $BOUNB. 
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Sample  Input  Decks 
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Several  sample  input  decks  are  presented  here  including  the  control 
cards  necessary  for  operation  on  the  CDC  6600  with  the  Scope  or 

NOS/BE  operating  system.  The  first  input  deck  (pp.  56  - 58)  is  the 
initial  input  stream  used  to  begin  the  three-dimensional  combustor 
simulation  presented  herein.  The  second  input  deck  (p.  59)  is  a restart 
of  a previous  run  from  restart  files  saved  on  magnetic  tape,  and  the  third 
input  deck  (p.  6o)  is  the  restart  deck  required  to  initiate  the  nitric 
oxide  (NDX)  chemistry  calctilation  frcm  a previous  run. 


I 

I 


I 
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Description  of  Printed  Output 

The  MINT  combustor  code  provides  detailed  printed  output  so  that  a 
user  can  monitor  the  solution  as  it  progresses  in  time,  and  examine  the 
flow  field  at  a user  specified  time- step- frequency.  The  initial  phase  of 
the  printed  output  (pp.  61-77)  contains  primarily  the  namelist  input 
described  above.  This  is  followed  by  the  initial  flow  field  print  of  which 
only  a sample  is  shown  on  pp.  78  - 8l,  and  then  a sample  periodic  time 
step  summary  is  shown  on  p.  82. 


The  first  line  of  output  is  the  heading  "MINT  COMBUSTION  CODE",  the 
second  line  will  be  either  "INITIAL  INRJT  STREAM"  for  a run  initiated  from 
cards  (IREST  = O)  or  "RESTART  OPTION  USED"  for  a run  which  is  a restart  of 
a previous  run  (iREST  >0),  The  next  line  of  output  appears  only  when  a 
new  restart  dump  is  to  be  written  (IREPUN  > O),  however  no  internal  position- 
ing of  restart  files  is  performed  by  the  code.  The  next  line  of 
output  reads  "DIMI-DIMIO  =",  followed  by  the  ten  reference  quantities: 

Re‘\  Pd/p^Uq’  Pd/°D^D’  (RePr)'-(U^h^), 

(ReSc)"^,  Bo’^,  and  aT^q^^,  respectively.  The  next  line  "T,  P,  H- TOTAL  =",  is 


followed  by  reference  values  for  these  quantities.  Then  namelist  $LIST2  is 
printed  by  subroutine  BVIV,  and  contains  the  grid  point  indices  shown  in 
Figs.  3 and  U.  If  IGEOM  < 0,  the  heading  "***^«geNERAL  GEOMETRY  COEFFI- 
CIENTS****" is  printed  followed  by  a tabulation  (for  all  Xj  and  Xj)  of  the 

metric  coefficients  and  their  derivatives  in  the  order  h2,  hp,  h2~^,  hT~^, 

Sh2/9n,  dhg/Ss,  B(h]_”^)/9n,  9(hj^"^)/9s.  The  three-point  central  difference 

operators  (6j^,  6^)  are  then  printed  for  directions  k = 1,2,3.  The  namelist 
$GRIDV,  which  contains  the  actual  coordinates  (X,  Y,  Z),  grid  spacing 
(DX,  DY,  DZ)‘,  and  the  Roberts'  grid  transformation  factors  (Xl,  X2,  Yl,  Y2, 
Zl,  Z2),  is  printed  by  subroutine  BVIV. 


The  remaining  namelist  variables  have  been  described  in  the  Input 
description  subsection.  The  namelist  $B0UND  is  printed  only  if  a new  case 
is  being  started.  Then  a single  line  is  printed  containing  the  nondimen sional 
time  steps  corresponding  to  the  explicit-stability  limits  due  to  viscous 
terms,  conduction  terms,  and  the  CFL  limit.  Subsequently,  the  namelists 
$READ1,  $READ2,  $READ3,  $READ4,  $READ5,  and  $CHEMD  are  printed  in  that  order 
followed  by  a labeled  "****  STARTING  FLOW  FIELD  ****"  print. 
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The  three  line  summary  output  px'inted  at  the  end  of  each  time  step  is 
largely  self-explanatory:  VISC,  STAB  and  CFL  indicate  the  ratio  of  the  actual 

time  step  DT  to  those  explicit  stability  limit  step  sizes;  SSTEST  is  the 
max~im\jm  relative  change  in  any  flow  variable  occurring  at  grid  point  IX,  LY, 
LZ;  RATIO  is  just  SSTEST/SSEPS;  DU,  IN,  DW,  DR,  DMF,  DTN,  DTKE,  DTED,  and  DNO 
are  the  maximum  relative  changes  in  the  variables  u,  v,  w,  p,  m^^,  T,  k,  e,  and 
“N0»  respectively.  This  is  followed  by  the  labeled  summary  o\itput  (NZ2  lines) 
from  subroutine  PRNTS. 

The  detailed  printout  obtained  at  a user  specified  time  step  increment 
is  clearly  labeled  as  follows:  U-VEL,  V-VEL,  W-VEL,  DENS  (p),  ENTH  (h  or  H), 
MFRACl  (m^),  MFEHI  ($)  if  chemistry  is  included,  TKE  (k),  TDISS  (e),  PRESS 
(p-PNOM),  TEMP,  VISC  (Mgff)*  DIV  (indicates  7-u),  ZMIX  (Z  = W^"^)  and 
CHMEX  (ch  = printer  plots  obtained  at  the  end  of  a run  are  also 

labeled  as  above;  TIME  PLOTS,  DIR-1  PLOTS,  DIR-2  PLOTS  and  DIR-3  PLOTS  are 
presented  in  that  order. 
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Figure  l(oontinued).  Flow  diagram  for  MINT  combustor  code 
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DROPLET  MASS  DISTRIBUTION  FUNCTION,  f(r) 


Table  I . References  quantities 


nominal  PRESSURE  PNOM  NONDIMENSONAL 


Table  II.  External  equation  and  variable  numbers 


GOVERNING  EQUATION 

lEQ 

EXTERNAL 

EQUATION 

NUMBER 

ASSOCIATED  DEPENDENT 
VARIABLE 

IV ARNO  UEQ) 
DEFAULT 
VARIABLE 
NUMBER 

ISOLVE  IIEQ! 
SOLUTION 
MODE  SENTINEL 

X,  - MOMENTUM 

1 

X,  - VELOCITY 
COMPONENT,  u 

1 

1 

Xj  - MOMENTUM 

2 

Xj- VELOCITY 
COMPONENT,  V 

2 

1 

X3  - MOMENTUM 

3 

X3  - VELOCITY 
COMPONENT,  w 

3 

1 

CONTINUITY 

4 

density,  p 

4 

1 

ENERGY 

5 

ENTHALPY,  H OR  h 

5 

1 

FUEL  SPECIES 
CONSERVATION 

6 

FUEL  MASS  FRACTION,  m, 

6 

1 

MIXTURE  FRACTION 
CONSERVATION 

7 

MIXTURE  FRACTION,  ♦’ 

7 

1 

NO  SPECIES 
CONSERVATION 

8 

NO  MASS  FRACTION, 

s 

3 

PARTICLE  conservation, 
CLASS  1 

9 

PARTICLE  FRACTION,  f , 

10 

2 

PARTICLE  CONSERVATION, 
CLASS  2 

10 

PARTICLE  FRACTION, 

11 

2 

PARTICLE  CONSERVATION. 
CLASS  3 

11 

PARTICLE  FRACTION.  13 

12 

2 

PARTICLE  CONSERVATION. 
CLASS  4 

12 

PARTICLE  fraction, 

13 

2 
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TABLE  II  (continued).  External  equation  and  variable  numbers 


GOVERNING  EQUATION 

lEQ 

EXTERNAL 

EQUATION 

NUMBER 

ASSOCIATED  DEPENDENT 
VARIABLE 

IVARNO  (IEQ) 
DEFAULT 
VARIABLE 
NUMBER 

ISOLVE  (lEQI 
SOLUTION 

MODE  sentinel 

X,  - RADIATION  FLUX 

13 

X,  - RADIATION  FLUX 
FUNCTION,  G, 

14 

2 

Xj-  RADIATION  FLUX 

14 

Xj  - RADIATION  FLUX 
FUNCTION.  Gj 

15 

2 

X3  - radiation  flux 

15 

X3  - RADIATION  FLUX 
FUNCTION,  G3 

16 

2 

TURBULENCE  KINETIC 
ENERGY 

16 

TURBULENCE  KINETIC 
ENERGY,  k 

8 

1 

turbulence  kinetic 

ENERGY  dissipation 
RATE 

17 

DISSIPATION  RaTE,€ 

9 

1 

Table  IV.  Default  values  of  lEQBC  (IBC,  lEQ) 


IBC 

1 

2 

3 

4 

5 

6 

1 

1 

1 

2 

2 

5 

3 

2 

2 

4 

1 

1 

5 

3 

3 

2 

4 

2 

2 

5 

3 

4 

2 

4 

2 

2 

4 

4 

5 

2 

2 

2 

2 

7 

3 

6 

2 

2 

2 

2 

5 

3 

7 

2 

2 

2 

2 

5 

3 

8 

2 

2 

2 

2 

5 

3 

9 

2 

2 

2 

2 

5 

3 

10 

2 

2 

2 

2 

5 

3 

11 

2 

2 

2 

2 

5 

3 

12 

2 

2 

2 

2 

5 

3 

13 

2 

4 

0 

0 

0 

0 

14 

0 

0 

2 

2 

0 

0 

15 

0 

0 

0 

0 

4 

4 

16 

2 

2 

2 

2 

5 

3 

17 

2 

4 

2 

2 

5 

3 

Table  VIII.  Default  variable  numbers  for  namelist  input 


IVAR 

EXTERNAL  EQUATION  NUMBER 
lEQINP  (IVAR) 

ASSOCIATED  INPUT  VARIABLE 

1 

1 

X,  - VELOCITY  COMPONENT 

2 

2 

Xj  - VELOCITY  COMPONENT 

3 

3 

X3  - VELOCITY  COMPONENT 

4 

4 

STATIC  PRESSURE 

5 

5 

TEMPERATURE 

6 

6 

FUEL  MASS  FRACTION 

7 

NONE 

TURBULENT  VISCOSITY 

8-11 

9-12 

PARTICLE  FRACTION.  CLASSES  1-4 

12 

13 

X,  - RADIATION  FLUX  COMPONENT 

13 

14 

Xj  - RADIATION  FLUX  COMPONENT 

14 

15 

X3  - RADIATION  FLUX  COMPONENT 
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